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Abstract 

In this talk, we briefly comment on Sweeny and Gliozzi methods, clus- 
ter Monte Carlo method, and recent transition matrix Monte Carlo for 
Potts models. We mostly concentrate on a new algorithm known as 
'binary tree summation'. Some of the most interesting features of this 
method will be highlighted - such as simulating fractional number of Potts 
states, as well as offering the partition function and thermodynamic quan- 
tities as functions of temperature in a single run. 

1 Introduction 

The Potts model |j is not only fascinating in relation to phase transitions but 
also an excellent model for constructing new algorithms for Monte Carlo simu- 
lation. Before the advent of cluster algorithms, Sweeny g] in 1983 already had 
an algorithm that has much reduced critical slowing down. Sweeny algorithm is 
a heat-bath algorithm in the Fortuin-Kasteleyn representation j| of the Potts 
model as a percolation problem for any positive real value of the number of Potts 
state q. Any thermodynamic averages in the Potts model representation can 
be translated into a percolation bond representation Q] . Sweeny algorithm can 
be efficiently implemented only in two dimensions. In three and higher dimen- 
sions, it becomes too expensive to simulate, because each move requires a step 
to identify if two sites belong to the same or different clusters. Nevertheless, it 
is perhaps the only algorithm that works for any real values of q. 

Very recently, Gliozzi |5| introduced a set of different transition probabilities 
than that of Sweeny's. In Gliozzi's rate, those links that are occupied by bonds 
or pair of unoccupied sites that connect the same cluster are resampled with 
probability p for occupation, and 1 — p for empty. For unoccupied sites that 
belong to two different clusters, the probability for occupation by a bond is p/q 
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and empty 1 — p/q. The resampling of a large fraction of the links appears 
to decorrelate the configuration efficiently. Gliozzi claims || that his algorithm 
does not suffer from critical slowing down. We have shown by explicit calculation 
of the correlation times that this is not true In two dimensions, we have 
good logarithmic divergence of the correlation time with system linear sizes 
for the Ising model. For two-dimensional three-state Potts model, we found 
dynamical critical exponent z w 0.5, and for the three-dimensional Ising model, 
the exponent is about 0.4. 

The Swendsen-Wang M and Wolff [0 algorithms work only for integer val- 
ues of q. They can be implemented very efficiently. For the Swendsen-Wang 
algorithm, each sweep takes CPU time of O(N), and for Wolff algorithm, it is 
proportional to the size of the cluster that is being flipped. They reduce but 
not completely eliminate critical slowing down at the critical temperature. 

The traditional methods mentioned above simulate at one fixed temperature 
T, thus the results are in terms of discrete points. A curve is obtained by 
connecting these points. The transition matrix Monte Carlo |J and other related 
algorithms [[l0| O, [ill [f3) generate continuous function of temperature or other 
variable from a single simulation, thus greatly enhanced efficiency. The multi- 
histogram sampling of Weigel et al [jl4| in the context of g-state Potts model 
should be mentioned here. The multicanonical simulation method of Berg jll], 
|f5[ also has the additional advantage of much reduced correlation times for 
first-order phase transitions. 

It is very attractive to have algorithms that do not have correlation at all. 
In fact, there are such algorithms, such as that of random percolation. Since 
the samples are drawn at random for each configuration, each configuration is 
independent of the previous configurations. The simple sampling method of Hu 
[ pi) has this property. Newman and Ziff |l7j gave an efficient implement that 
work by exact weighting of the probability p so that after the simulation, the 
whole curve Q of some physical observable versus p can be computed. In this 
paper, we present some attempts [|f8| that generalize Newman-Ziff method to 
general Potts model of q states. These algorithms all have the common feature 
that bonds are added to the system one by one. Once there are on the lattice, 
they do not move. It is a noncquilibrium process in the sense that there is no 
(Monte Carlo) time translation invariance for the process. The rest of the paper 
is organized as follows: in the next section, we briefly review the Newman-Ziff 
method, followed by generalization of the partition function to Potts model. We 
then discuss a number of algorithms, and show their efficiencies and also point 
out their shortcomings. We conclude in the last section. 

2 Newman-Ziff-type methods 

Consider the simple bond percolation problem. We generate configurations by 
putting bonds at every nearest neighbor link of a (/-dimensional hypercubic 
lattice with probability p, and absence of bond with probability I — p. Any 
average of physical quantity Q(r) of configuration T can be computed exactly 
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from 

(Q)=^Q(lV(l-p) M - b , (1) 
r 

where the summation is over all the 2 M configurations; b is number of bonds 
present, and M = Nd is maximum possible number of bonds. 

In a standard simulation, one fixes a value of p, and generates configura- 
tions by visiting each link and putting a bond with probability p. Estimates 
of (Q) are obtained by taking sample means of quantity Q(T). Newman and 
Ziff |nj instead considered a computation of Qb for each value of b, and then 
reconstructed from Qb the function (Q) of p. From Eq. (|IJ), we have 

(Q) = I2Y.Q(r b )p\i-pr- b = fy ( i q 6> (2) 

6=0 T b 6=0 '' 

where Qb is defined as average over all the M\/(b\(M — £>)!) configurations Tb of 
exactly b number of bonds, 

Qi= bl M_^ EQ(ri) (3) 

Newman and Ziff proposed the following algorithm to evaluate Qb- Each 
sweep of the lattice starts with empty lattice, putting bonds one by one at 
random over these that is empty. This will generate M + 1 configurations with 
0, 1, • • •, b, ■ ■ •, M number of bonds. From them, we calculate Qb for each b. 
With the help of Hoshen-Kopelman algorithm jL9) , the whole calculation of one 
sweep can be done in CPU time of O(N). Although configurations within a 
sweep are highly correlated, each sweep is an independent one. Very accurate 
percolation threshold was obtained this way. 

To generalize this procedure to Potts model, we note that the analogous 
equation to Eq. (||) is 

M M 

(Q) = E £ Q^il - p) M ~ b q N ^ = $> l (l - p) M - b C b Qb, (4) 
6=0 T b 6=0 

where N C (T) is the number of clusters of the configuration T, and 

£^ c(r6) ' Q b = -£Q(r6)^. (5) 
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The probability p is related to temperature T by p = 1 — exp(— J/(fcT)). When 
q ^ 1, Cb is no longer known exactly, we must compute by Monte Carlo sim- 
ulation; the quantity Qb is now an average over the configurations of a given 
number of bonds b distributed not uniformly but according to q Nc ( Fb \ This 
makes the simulation much harder. 
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3 Algorithms for computing q, and 



3.1 simple surviving and dying process 

The following algorithm, although very inefficient, generates correct probability 
distribution for the samples, and introduces the way c& and in general any 
can be computed. 

Starting with an empty lattice, one sweep consists of repeated application 
of the following steps until the process dies: 

1. Pick an unoccupied neighbor pair at random for the next bond. 

2. If inserting a bond 

(a) does not change the cluster number (AiV c = 0), accept the configu- 
ration; 

(b) merge two clusters, so that the cluster number decreases by 1, (AN C = 
— 1), accept the configuration with probability l/q, or reject the con- 
figuration and terminate the process (and begin the next sweep from 
an empty lattice). 

3. Take statistics of the survival configurations (with equal weights). 

The probability distribution of the process at b number of bonds is proportional 
to q Nc ( r \ The configurations that have equal number of clusters from b to 
6+1 bonds have the same probability, but these that merge clusters appear 
with a probability smaller by a factor l/q. The cumulative effect gives the 
desired probability distribution. Unfortunately, since the number of samples is 
exponentially small for large b, the algorithm is of only conceptual use. 

The number Cb can be related to sample average of the conditional survival 
probability. Let uq be the number of empty links that connect same cluster, 
and ni be the number of empty links that are on different clusters. Then the 
conditional survival probability is fis}| 



To speed up the simulation, we prevent the process from dying, but the price we 
have to pay is that we need to give weight to samples. Specifically, we follow the 
method of ./V-fold way: we choose a type-0 class with probability uq/ (no -\-n\/q) 
or type-1 class with probability n\/qj(nQ + ri\jq). Once we have decided the 
class, we pick an empty link from among no bonds (for type-0 class) or among 
ri\ bonds (for type-1 class) at random. The probability distribution of such a 
process is not exactly what we wanted, but rather it is 





3.2 N-fold way 



-Ppath(rb) 



1 



,N c (r b )-N a (r ) 



(7) 
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where A(T) = no(T) + m(T)/q, and N C (T ) = N is the number of clusters 
of the empty lattice. To compute the desired average we must weight with 

J2« Nc{rb) Q( T *) E 4 ot Q(r 6 ) 

. . T;, MC samples 

Tf, MC samples 

Due to fluctuations in the weights, inefficiencies are unavoidable. In fact, the 
above method can give realizable results only for small sizes of order 10. 

3.3 binary tree summation 

To overcome the above problem, we propose the following method |lq] named 
binary tree summation Monte Carlo. For each sweep, the algorithm consists of 
two parts: the simulation part and summation part. In the simulation part, 
we always pick a type-1 link (unoccupied link that bridges two clusters) at 
random, thus always merge clusters in each step. The probability distribution 
of the configuration sequence generated is 

p(r -> r t ->■••-> r<) = 1 1 (9) 

m(r )ni(ri) ni(rvi) 

In the summation part, we consider all possible ways of inserting type-0 bonds 
in between each merge step of simulated configurations. The configurations are 
weighted in such a way to realize the desired probability distribution, P(T) oc 
qN c (r)_ Specifically, the weight of a configuration which is specified by number 
of bonds b and sequence number in the simulation i, is the product of factor 
ni/q (for each merge step) and no (for each insertion of type-0 bond that does 
not merge clusters). This weight w(b,i) can be computed recursively for all b 
and i with computer time of 0(N 2 ), by 

w(b+l,i) = w(b,i)n (b,i) + w(b,i - l)n x {b,i- l)/q, (10) 

where ni(b,i) = ni(i), rio(b,i) = rio(i) — b + i, with the starting condition 
w(0,i) = Sifl and the constraint no(b, i) > 0. The value w(b,i) is nonzero 
only for b > i. The computation of the weights w(b,i) is similar to a recursive 
computation of the binomial coefficients. The final statistics of a quantity Q 
can be computed as 

N-l 



Q b = (w b )- i {Y,™(b,im)), (ii) 



i=0 



where W b — X^o 1 w (bii)> Q(i) 1S the quantity at the i-th cluster merge step. 
The ratio (b + l)c b +i / ((M — b)c b ) can be computed from the expectation value 
of (n + ni /q)/(M-b). 
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Figure 1: Specific heat for the two-dimensional Ising model computed by binary 
tree summation method. The insert gives the relative error with respect to exact 
result. The system size is 32 x 32 with 10 6 Monte Carlo sweeps. 



Note that the number of states q enters into the simulate only through the 
weights and is not used in sampling, we can use any real or complex value for 
q. Moreover, if we collect appropriate histogram (of 0(N 2 ) entries), we can 
reconstruct the equilibrium average for any p and q through a single simulation. 
The details of the algorithm and some applications will be discussed elsewhere 

In Fig. |1|, we present a plot of specific heat of the Ising model on a 32 x 32 
lattice. Errors in comparison with exact results are indicated in the insert. 
For this lattice size, the errors are quite small for all values of T. In Fig. ||, we 
plot the average relative errors in q,, defined by 

(12) 

The quantity c& plays the role of density of states as considered in refs. p2l 
and p3[ . For the binary tree summation data, we have rescaled the error by 
a factor y/t cpu /1.9 so that the comparison of errors is on an equal footing of 
given amount of CPU time. As we can see from the plot, the errors for small 
lattices are much smaller than that of both the multi-canonical simulation and 
transition matrix Monte Carlo method 0. However, as the system size becomes 
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Figure 2: Average relative errors in c& after 10 6 sweeps for the two-dimensional 
Ising model computed by binary tree summation method (circles), and errors 
in density of state n(E) from multicanonical method (squares) and a two-stage 
transition matrix Monte Carlo method (triangles). The data are from ref. |L8j 
andg. 

larger than 32 x 32, the error becomes comparable and less favorable to the other 
methods. 

The decreased performance of binary tree summation method for large sys- 
tems is partly due to the 0(N 2 ) nature of the algorithm, and perhaps more 
importantly, it is not an importance sampling method. The fluctuation of the 
total weights Wb m binary tree summation method is much smaller comparing 
to the iV-fold way. Nevertheless, the weights introduce extra errors and become 
important for large systems. The q = 1 percolation problem do not have such 
problem as in that case the total weight for each b is a constant. A "fudge 
factor" may be introduced in the probability of choosing a type-1 bond with the 
aim to reduce the total weight fluctuation. We are still investigating on such 
possibilities. 

4 Conclusion 

We first compared the Sweeny and Gliozzi rates for the simulation of Potts 
models. The Sweeny and Gliozzi rates give the same dynamics and also do not 
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completely eliminate critical slowing down. We then introduced Newman and 
Ziff method for simulating percolation problem. Our central goal is to simulate 
Potts models in the Fortuin-Kasteleyn representation in the spirit of Newman 
and Ziff. Binary tree summation method is our best attempt in this direction. 
Although the samples are generated without correlation between sweeps, we 
lost the property of importance sampling. In this sense, we have not completely 
solved the problem, thus the methods are not efficient for large systems. We 
hope the present work will stimulate research for algorithms that have zero 
correlation, yet realize importance sampling. 
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